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Abstract 
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long-range tail in the momentum spectrum. These are accounted for in detail by discussing the natural 
orbitals and their occupations. Our discussion is complemented by an analysis of the two-body correlation 
function. 
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I. INTRODUCTION 



For more than a decade, trapped ultracold atoms have been a vastly expanding field of research 
|[liy,y,y]- Their low thermal fluctuations, and the fact that one can virtually design both the ex- 
ternal as well as the inter-particle forces via electromagnetic fields, make them an ideal simulation 
tool for various phenomena, ranging from condensed-matter physics to quantum information. A 
particularly intriguing aspect has been the system's dimensionality. For example, if its transverse 
degrees of freedom are frozen out such that an effective one-dimensional description becomes 
possible, then the effective interaction strength of the system may be tuned at will from a weakly 
correlated to a strongly interacting one by merely changing the lengthscale of the transverse con- 
finement lO]. In the case of infinite repulsion — the so-called Tonks-Girardeau gas — the system 
may even be mapped to that of an ideal Fermi gas [j^ and exhibits striking features, including the 
reduction of off-diagonal long-range order vj] or a very distinctive momentum spectrum Isl. This 
has also been verified experimentally IgiflOl. 

However, the Tonks-Girardeau limit requires low densities and is thus amenable only for few 
atoms, typically ~ 10 — 100. Besides that, there is the factor of computational accessibility — 
which includes the fact that a small number of atoms greatly facilitates the intuitive understanding 
of the mechanisms underlying the physics of large systems. But there is yet another argument for 
considering few-atom systems: They permit a much higher level of control. There is no thermal 
cloud, as for large N, associated with decoherence and energetically dense excitations, but a pure 
quantum system. 

However, the same reasons that make few atoms so interesting also render their solution com- 
putationally cumbersome. There have been some detailed analyses of the Tonks-Girardeau limit in 
a harmonic trap which is greatly simplified since the solution is semi-analytic. Needless 

to say, the complementary borderline case of weakly interacting bosons is also well understood: in 
this limit, the Gross-Pitaevskii equation becomes valid, which assumes that all particles condense 
into a single one-particle orbital (see, e.g., [1]). Comparatively little is known about the transition 
between these extremes, though. Many key features are nicely illustrated on the analytic solution 
of two bosons in a harmonic trap lisl or other simple models 1 14, jj]. Based on a multi-orbital 



mean-field approach valid for arbitrary traps, it has recently been demonstrated that there exists an 
intriguing pathway for the above transition Ilia] . While that ansatz captures the essential features 
of that evolution, it is desirable to investigate it from a rigorous many-body perspective, but still 
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irrespective of the trap geometry. Most of the studies so far have relied on the assumption of only 
few contributing orbitals IitI [isl . while a recent exact diagonalization has focused on a simple 
harmonic trap lll9ll . For the case of a double- well trap, the authors have previously investigated the 
ground state of few interacting bosons in the entire regime between the Gross-Pitaevskii limit and 
fermionization for different well depths, but chiefly from the perspective of the density profiles 

Q. 

The aim of this paper now is to address the role of correlations in these systems, in particular 
the relation of the fragmentation mechanism to the concept of long-range order and the momen- 
tum density. This way we not only bridge the gap between the two borderline regimes of zero 
and infinite repulsion but also extend our previous stud y. O ur approach is based on the Multi- 
Configuration Time-Dependent Hartree method which we use to study the numeri- 
cally exact ground state of few bosons. 

This article is organized as follows. In SecHIl the model is introduced and some key quantities 
for the description of correlations are discussed. Section Hill contains a concise introduction to the 
computational method and how it can be applied to our problem. In the subsequent section the 
one-body correlation aspects are studied. More concretely, we present the full one-body density 
matrix in the context of long-range order in Sec. IIV Al which is subsequently explained in terms 
of its eigenvectors and their populations (IIV BD . Their connection to the momentum density will 
then be clarified in Sec. IIV CI We complement our investigation by going beyond the one-particle 
picture in looking at the two-body correlations in Sec. |Vl which is rounded off by relating our 
results to some common approximation schemes, the two-mode model |^] and the multi-orbital 
mean field. 



II. THEORETICAL BACKGROUND 
A. The model 

In this article we investigate a system of few interacting bosons (A^ = 2, . . . , 6) in an external 
trap. These particles, representing atoms with mass M, are taken to be one-dimensional (ID). 
More precisely, after integrating out the transverse degrees of freedom and rescaling we arrive at 
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the model Hamiltonian (see ll20ll for details) 




i i<j 



where h = + U (x) is the one-particle Hamiltonian with a trapping potential U, while V is the 
two-particle interaction potential with low-energy scattering length oq, taken to be the effective 
interaction ol 



Here a harmonic transverse trap potential with oscillator length a±_ was assumed. Moreover, the 
well-known numerical difficulties due to the spurious short-range behavior of the delta-function 
potential 5{x) are alleviated by mollifying it with the normalized Gaussian 



which tends to 5 as (j ^ in the distribution sense. We choose a fixed value a = .05 as a trade-off 
between smoothness and a short range. 

B. Fragmentation: key aspects 

Although our approach equips us with the full solution of the system — here, the ground-state 
wave function — this solution obviously still needs to be related to the concrete physical questions. 
Penrose and Onsager suggested a criterion connected with the one-body density matrix, which 
will be laid in what follows. 

As is well-known, the knowledge of the wave function ^ is equivalent to that of the density 
matrix = |^)(^|. To the extent that we study at most two-body correlations, it already suffices 
to consider the two-particle density operator 



whose diagonal kernel p2{xi, X2) gives the probability density for finding one particle located at 
Xi and any second one at X2- For any 1 -particle operator, of course, it would be enough to know 
the one-particle density matrix pi = tr2P25 so that the exact ground-state energy may be written as 
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P2 = tr3..jv|^)(^ 



(1) 



E = NtT{pih) + 



N{N-1) 
2 



tr(p2^). 
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Consider the spectral decomposition of the one-particle density matrix 



Pl = ^na\(f)a){(j), 



(2) 



a 



where Ua E [0, 1] is said to be the population of the natural orbital (pa- If all n'^ = UaN G N 
(Xla ~ density may be mapped to the (non-interacting) number state |nQ, ri'^, . . . ) 

based on the one-particle basis {0a}; for non-integer values it extends that concept. In particular, 
the highest such occupation, no, r nay serve as a measure of non-fragmentation, a criterion put 
forward by Penrose and Onsager llisll . For no = 1, a simple condensate is recovered. This is 
the well-known borderline case of the Gross-Pitaevskii eq.: as (7 ^ 0, pi ^ |0o)(0o| L26] and 
P2 = Pi ® Pi, so that the interaction above can be replaced by a mean field V = tr(piV^). 

The 1 -particle density matrix may be viewed not only from the perspective of its spectral de- 
composition, but also in terms of its integral kernel pi{x,x') = {x\pi\x') = pi{x',x)*. Since 
the density matrix is non-negative, so is the 1-particle density p{x) = pi{x,x). As opposed to 
that, the off-diagonal part will be complex in general (though in this paper a real representation is 
employed). It is therefore certainly not an observable in its own right. Nonetheless, it is of some 
interest as it gives us access to all one-particle quantities, also non-local ones such as the mo- 
mentum density p(A;) = 27r(A;|pi|A;) = ^^na|0a(A;)p, which can be related to the density matrix 



It is reflection symmetric if pi is real symmetric. Moreover, it can be understood as the Fourier 
transform of the integrated 'off-diagonal' correlation function 111] 



with 7(r) := J dRpi{R + — |). Note that 7 is again generally complex and reflection 
symmetric, while 7(0) = 1. From this, it becomes clear that the off-diagonal behavior encoded in 
7 has a 1-1 correspondence to the momentum distribution. More specifically, the short-distance 
behavior determines the high-/c asymptotics, which for a delta-type interaction V(x) = g5{x) in 
the limit g ^ 00 has been shown to display the universal decay p{k) = 0{k^^) [8]. Conversely, 
the off-diagonal asymptotics r ^ 00 relates to the low-A; regime. This, however, depends on the 
nature of the external potential. For a translationally invariant system, it has been argued that Bose 
condensation were, equivalentto off-diagonal long-range order, i.e. 7(r) = 0(1) aim . In the same 



via 
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context, but in the limit (? cxo, it has in turn been shown that 7(r) = 0(r^^/^), which implies an 
infrared divergence p{k) ~ c/ Vk as A; ^ 1^]. 

The above limit g ^ oois commonly referred to as the Tonks -Girardeau limit of ID hard-core 
bosons, or also as their fermionization. This lingo finds its justification in the Bose-Fermi map 
that establishes an isomorphy between the exact bosonic wave function \1'+ and that of a 
(spin-polarized) non-interacting /erm/on/c solution \1'q , 

where A{Q) = ]\i<:j sgn(xj — Xj) and Q = (xi, . . . , xn)'^- In particular, the ground state reduces 
simply to the absolute value of \I/q , which makes it tempting to think of the hard-core interaction 
— i> oo as mimicking the exclusion principle. As the free fermionic solution is easily accessible, 
this theorem has proven very fruitful in a wide range of applications. For our purposes, the most 
relevant one is the solution of bosons in a harmonic trap llli 

VI/+ (g) oc e-l«l'/^ n l^^-^^-l' 
i<«<i<-'v 

which illustrates the characteristic short-distance correlations. 



III. COMPUTATIONAL METHOD 



Our goal is to investigate the ground state of the system introduced in Sec. HU for all relevant 
interaction strengths in a numerically exact, i.e., controllable fashion. This is a highly challeng- 
ing and time-consumi ng task, and only few such studies on ultracold atoms exist even for model 



systems (see, e.g., u2 
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1911'). Our approach relies on the Multi-Configuration Time-Dependent 
Hartree (MCTDH) method 921 . Iisl lioi . primarily a wave-packet dynamics tool known for its out- 
standing efficiency in high-dimensional applications. To be self-contained, we will provide a 
concise introduction to this method and how it can be adapted to our purposes. 

The underlying idea of MCTDH is to solve the time-dependent Schrodinger equation 

vi/(g,o) = vi/o(g) 

as an initial- value problem by expansion in terms of direct (or Hartree) products $ j: 



(3) 



"1 



/ 



(4) 



using a convenient multi-index notation for the configurations, J = (ji . . . jf), where / denotes the 
number of degrees of freedom and Q = (xi, . . . ,Xf)'^. The (unknown) single-particle functions 
(p^j^J are in turn represented in a fixed, primitive basis implemented on a grid. For indistinguishable 
particles as in our case, the sets of single-particle functions for each degree k = 1, . . . , are of 
course identical (i.e., we have ipj^, with < n). 

Note that in the above expansion, not only the coefficients Aj are time-dependent, but so are 
the Hartree products $j. Using the Dirac-Frenkel variational principle, one can derive equations 
of motion for both Aj, $j ll23ll . Integrating this differential-equation system allows one to obtain 
the time evolution of the system via Let us emphasize that the conceptual complication above 
offers an enormous advantage: the basis t)} is variationally optimal at each time t. Thus it 

can be kept fairly small, rendering the procedure very efficient. 

It goes without saying that the basis vectors <l>j are not permutation symmetric, as would be 
an obvious demand when dealing with bosons. However, the symmetry can be enforced on \l/ by 
symmetrizing the coefficients Aj, even though this turns out to be unnecessary as long as identical 
single-particle functions are employed. 

The Heidelberg mctdh package which we use, incorporates a significant extension to 



the basic concept outlined so far. The so-called relaxation method 113 Ih provides a way to not only 
propagate a wave packet, but also to obtain the lowest eigenstates of the system. The underlying 
idea is to propagate some wave function \E'o by the non-unitary e"^"" {propagation in imaginary 
time.) As r ^ oo, this automatically damps out any contribution but that stemming from the true 
ground state | J = 0), 

e-^-vl/o = ^e-^-'^|J)(J|M/o). 

In practice, one relies on a more sophisticated scheme termed improved relaxation. Here 
{'^\H — E\^!) is minimized with respect to both the coefficients Aj and the configurations $j. 
The equations of motion thus obtained are then solved iteratively by first solving for Aj (by diago- 
nalization of {{^j\H\^k)) with fixed $j) and then propagating <l>j in imaginary time over a short 
period. That cycle will then be repeated. The improved-relaxation method is outlined in Ref. Illi : 



a more comprehensive account is also available ||3^ 

As it stands, the effort of this method scales exponentially with the number of degrees of free- 
dom, . This restricts our analysis in the current setup to about = 0(10), depending on how 
decisive correlation effects are. If these are indeed essential, it has been demonstrated Kflll that at 
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Figure 1: Sketch of the model potential U{x) = + h5w{x), consisting of a harmonic trap plus a 
normalized Gaussian of width w = 0.5 and barrier strengths h = 2,5, 10. 

least n = N orbitals are needed for qualitative convergence alone, while the true behavior may 
necessitate many more. As an illustration, we consider systems with ~ 5 and use n ~ 15 
orbitals. By contrast, the dependence on the primitive basis, and thus on the grid points, is not as 
severe. In our case, the grid spacing should of course be small enough to sample the interaction 
potential. We consider a discrete variable representation l33i of 95 harmonic-oscillator functions, 
which is equivalent to 95 grid points. 



IV. ONE-BODY CORRELATIONS 

As in Ref. Ilfll . we consider the ground-state properties of bosons in a double- well trap modeled 

by 

U{x) = -x^ + hSw{x). 

This potential is a superposition of a harmonic oscillator (HO), which it equals asymptotically, 
and a central barrier which splits the trap into two fragments (Fig. [Hi. The barrier is shaped as a 
normalized Gaussian 5w of width w = 0.5 and 'barrier strength' h. 

For h = 0, the case of interacting bosons in a harmonic trap is reproduced. In Ref. l2fli . we have 
described the transition from a simple, weakly interacting condensate (g 0) to fragmentation 
and finally the Tonks -Girardeau limit (g oo). As /i ^ oo, the energy barrier will greatly exceed 
the energy available to the atoms, and we end up with two isolated wells. A larger g then affects 
only the fragmentation within each of these wells. In between, there is an interesting interplay 
between the 'static' barrier (h) and 'dynamical barriers' in the form of inter-particle forces (g). 
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Figure 2: 1-particle density matrix pi{x,x') for = 5 bosons. Top row: harmonic trap, bottom row: 
double well (barrier height h = 5). Results are shown for the interaction strengths g = 0.4, 4.7, 194 from 
left to right. 

which has been analyzed mainly from the viewpoint of spatial densities. 

We now seek to extend that investigation to non-local properties, such as the off-diagonal den- 
sity matrix (Sec. lIV Ab and the momentum density dlVCI) . so as to attain a deeper insight into the 
nature of the above transition. The role played by correlations will also be highlighted by showing 
how all the natural populations evolve as a function of g. 



A. One-particle density matrix and long-range order 

The 1-particle density matrix pi contains all the information about the one-particle aspects of 
the system, and serves as a good measure for the degree of fragmentation. In this section, we will 
analyze it from the most immediate perspective, i.e., we investigate its integral kernel pi(x, x'). 
Although it is not an observable in itself, being generally complex- valued, it is indirectly accessible 
e.g. via interferometry experiments yj]. More importantly, the density matrix relates to salient 
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questions such as off-diagonal long-range order or the momentum distribution (see 111 BD . Insofar 
it serves as a good starting point for our discussion of fragmentation in double- well systems. 

In FigureEl the fragmentation transition as reflected in pi (x, x') is visualized for = 5 bosons 
in a harmonic trap {h = 0, top row) and a double well of barrier strength h = 5 (bottom). In the 
harmonic case, the system starts at g = with a direct-product state -0 = 0o i-C-, with a density 
matrix pi(x, x') = (j)o{x)(j)Q{x') oc e~^^e~'^^^'^ in terms of r = x — x' and 2R = x + x'. From this 
point of view, the system does not exhibit off-diagonal long-range order, which is simply rooted in 
the fact that it is spatially bounded. Of course, it is nonetheless in a coherent state and thus features 
weak long-range order in that — x) ~ a/ p{x)p{—x) as x ^ oo. This property persists so 
long as the correlations induced by the interactions are weak enough for the system to remain 
in such a single-particle state (the Gross-Pitaevskii regime), such as for g = 0.4. For g = 4.7, 
however, the symmetry in R and r breaks up. The density profile p(x) = Pi{x, x) flattens, and one 
can see that the off-diagonal range is somewhat extended, too. However, as g is increased further, 
the support of pi(x, x') will concentrate more and more in the central region {x = x'}, where the 
typical fermionized profile is recovered (cf. g = 194). By contrast, the off-diagonal contributions 
will be washed out, indicating the decoherence of the system. Still it is noteworthy that even in 
this limit, a rest of coherence is preserved in a faint checkerboard pattern. 

For the double well (h = 5; bottom row), the situation is apparently different. As always, the 
system exhibits coherence to begin with (cf. g = 0.4), only that the orbital is now delocalized in 
both minima ±xo and may be written as (po{x) = c[ip{x — xq) + (fix + xq)]. Unlike the harmonic 
case, the off-diagonal range is not initially increased but directly destroyed upon switching on g. 
While for g = 4.7, the density matrix pi (x, x') may still be thought of as pertaining to two separate 
subsystems, it eventually reaches the Tonks-Girardeau limit (g = 194), where the only essential 
difference toward h = consists in the density suppression at x, x' = 0. 

B. Natural orbitals and their populations 

While, in principle, the full density matrix pi(x, x') as studied in the previous section contains 
all the information about fragmentation at the one-particle level, it is somewhat less amenable to 
intuition. A handier criterion is offered by its spectral decomposition Q in terms of its natural 
orbitals (f)a and their populations n„, telling us how close the system is to a pure one-orbital state. 
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Figure 3: Natural populations na{g) (a < 13) for iV = 4 bosons in a harmonic trap (a) and in a double well 
with barrier height /i = 5 (b) , /i = 10 (c). 



1. Natural populations as a measure of fragmentation 

Figures Ha-c) show typical plots of the natural populations as the interaction is increased, 
{na{g)}, for four bosons and h E {0, 5, 10}. Starting from uq = 1 for the non-interacting case, the 
lower lines rise steeply until they end up saturating in a fermionized state at (7 — > oo. Note that this 
pattern is to some degree universal. In other words, it is roughly detached from the specific shape 
of the trap, i.e, from what the underlying orbitals look like. This indicates why the set {ua} lends 
itself as a handy criterion for fragmentation. The details of the system are essentially encoded in 
(i) the exact sequence of Ua in the Tonks-Girardeau limit, and (ii) in the transition between the two 
extreme regimes = and g ^ 00. 

For the harmonic oscillator (h = 0), the plot reveals a relatively simple hierarchy. The value 
of no decreases smoothly to its Tonks-Girardeau limit A^~°-^^ llli . All the remaining populations 
increase dramatically up until (7 ~ 10, and accumulate in a more or less equidistant spacing (on a 
log scale). But even the next-to-dominant weight rii is nowhere near the 'condensate' fraction uq; 
the obvious gap between these two reflects the difficulty to observe fragmentation in the harmonic 
oscillator as compared to h > 0. Note that the group of lines {uq, . . . , nN_i} reveals a discernible 
separation from the lines below. This clearly relates to the finding that the convergence of the 
energy E(n) as a function of the number of orbitals, just to give an example, gets strikingly better 
once n It is the accumulation of points na{g) that makes for the utter slowness of 

true convergence pointed out in that work. 
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For a barrier with height /i = 5, a little more structure can be identified in the line sequence 
na{g). The accumulation persists, but at least the more populated orbitals a seem to come in 
groups of two. This will become clearer when looking into the natural orbitals. Even more strik- 
ing is the behavior of the second orbital's population, rii. It increases with g much more rapidly 
than all others, and it becomes comparable with uq already for modest ~ 5. This scale separa- 
tion between the pair no/i and the rest is in sharp contrast to the HO case. It gives a qualitative 
justification of the two-mode approximation widely used in double-well systems. To make these 
points even clearer, we have plotted the results for a much higher barrier, /i = 10. Here ni 'jumps' 
almost instantaneously {g <^ 1) to a value of order |, whereas the remaining occupations only 



catch up for (7 ~ 5. It is in that regime that the 2-mode model works brilliantly (see also Sec. lVBI) . 

The reason why fragmentation is facilitated when the central barrier is raised is intuitively clear. 
The particles' tendency to separate due to repulsion is usually obstructed by the higher costs of 
kinetic and potential energy. The potential-energy barrier creates an additional incentive for the 
bosons to separate. This has also been argued on more quantitative grounds (see, e.g., |24fl'). In a 
naive single-particle picture, the energy gap A in a double well between anti- and symmetric state, 
= c[ip{x — xq) ± (^{x + xq)], vanishes as /i — > 00. It is thus far easier for the interaction to 
bridge that gap for larger barriers, in particular compared to the gap for /i = 0, A = 1. 

Our final remark concerns the dependence on the atom number A^. For odd A^, two features of 
the even-A^ picture will differ. First, the second mode is less relevant (e.g. for = 5, ni|g=4 7 = 
0.16). What is more, the separation between the first A^ populations and all others was found to 
be much smaller. This backs up the intuitive notion that, for odd A^, fragmentation is seemingly 
impeded [20]. 



2. Natural orbitals 

Even though the natural orbitals (0a) are not of direct physical importance, they are a valuable 
tool to gain some insight into the process of fragmentation, as they determine both the spatial 
density matrix pi (x, x') as well as the momentum density p, to be discussed in the following sub- 
section. In the uncorrected case (7 = 0, the system is in a number state |A^, 0, . . . ) and thus the 
natural orbitals coincide with the single-particle eigenstates. Since V is a continuous perturbation, 
the orbitals 0^ will be somewhat distorted in the course of increasing g. For small enough g — i.e., 
in the Gross-Pitaevskii regime — that modified 0o will suffice for an accurate description. Con- 
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Figure 4: Natural orbitals for the case = 4. Top row: 'condensate' orbital 4>q for /i = (a), h = 5 (b). 
Bottom: 01 for /i = (c), h = 5 (d). Each subfigure shows plots for the representative interaction strengths 
g = 0.4 (weak repulsion), g = 4.7 (intermediate regime), and g = 7A (fermionization limit). 

versely, if correlations are sufficiently influential, many orbitals will contribute to pi, and studying 
their interplay will illuminate our results on the density matrix and the momentum distribution. 

Harmonic trap (h = 0) For the harmonic trap (Figs.|4^,c), the initial HO function 0o is only 
slightly flattened in the Gross-Pitaevskii regime (cf. g = 0.4). The onset of fragmentation not 
only smears out the lowest orbital, but also admixes an antisymmetric HO-type orbital 0i. In the 
fermionization limit, it is astonishing that already (po exhibits all the features of the fermionized 
density profile p(x), that is, pronounced humps mirroring the spatial isolation of the atoms. 
This is intelligible given that 0o still has a dominant weight, which ought to be contrasted with 
the philosophy of multi-orbital mean-field schemes (see, e.g., Ref. [16]), where that pattern is 
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produced by N spatially localized orbitals of equal population. That issue will be discussed in 
more detail in Sec. IV Bl 

Interesting as the orbitals may be in their own right, they also prove helpful in clarifying the 
decoherence found in Sec. lIV Al The onset of fragmentation, as for g = 4.7, leads to a broadened 
diagonal profile pi(x, x), but not equally so for the off-diagonal part. That is simply because the 
(pa have alternate parity (—1)'', and thus the admixture of another orbital leads to pi(x, —x) = 
J2a \4'a{x)f- Hcncc the fragmentation into different orbitals tends to deplete the off- 

diagonal as compared to the diagonal density. For g = 4.7, this effect is still tiny as ni ~ 0.1 
only, and therefore outweighed by the altogether extended support of 0o- However, as more and 
more orbitals are mixed, as is the case in the fermionization limit (see g = 74), this decoherence 
attains its full impact. We remark that the faint checkerboard pattern (Fig.|2l) is still rooted in the 
dominance of the lowest orbital, uq ^ A^"*^ "^^. 

Double well (h = 5) In the case of a central barrier (Figs. |lb,d), the natural orbitals in 
the non-interacting limit will again be the single-particle eigenstates, approximately the symmet- 
ric and anti-symmetric states (l)±. In the Gross-Pitaevskii regime (g = 0.4), the lowest orbital 
is only marginally flattened due to interactions, but a tiny reduction of the off-diagonal peaks 
Pi(±Xo, T^o) hints already at a minor admixture of the antisymmetric 0i. For g = 4.7, fragmen- 
tation has set in, not only smearing out the orbitals 0o/i — and thus the diagonal profile — but along 
the way washing out the off-diagonal long-range order almost completely. As emphasized before, 
the fermionization pattern tends to be generic for different h, which reflects both in the density 
matrix as well as in the natural orbitals. 



C. Momentum density 

The discussion so far focused on rather abstract aspects of the one-body correlations. Yet it can 
help us cast a light on an experimentally more amenable quantity, the momentum density. 
Jiarmonic trap (h = 0) For this case, the momentum distribution has recently been computed 



(13]; see also Ref. ll34l] ). yet we plot it in Figure |5] for comparison. It evolves from a Gaussian 
p{k)/2'K = TT^^/^e^'^^ at g = (with a maximum at p(0) = .35 . . . ) to a slightly sharper peak, 
here depicted for g = .4. This squares with the broadened natural orbital 00 in that regime, as 
found in Sec. lIVB'2l In the same light one sees that for g = 4.7, where fragmentation has set in, 
the peak at = is even more pronounced, while p(/c) has also developed a long-range tail. Both 
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Figure 5: Momentum density p{k) for (a-b) N = 5 bosons in a double-well trap of barrier height h. (left: 
h = 0, right: /i = 5); shown are the interaction strengths g = 0.4, 4.7, 15. (c) The case N = A, h = 10 for 
g = 0.01, 0.2, 4.7, 74. 

observations are easily accounted for. The k = behavior, for one thing, was argued to correspond 
to the off-diagonal long-range behavior of pi{x, x') in Sec. Ill Bl This fits in with our observation 
that the off-diagonal range was indeed extended in that ^f-regime, as seen in Fig.|2l 

The asymptotics A; ^ oo is in turn determined by the short-range interaction, which is known to 
culminate in the tail for g oo. This latter consequence is in fact confirmed here (see g = 15). 
Moreover, notice that the A; = peak is bound to diminish. In other words, the momentum 
spectrum is redistributed toward higher k, in accordance with the reduction of off-diagonal long- 
range order. This fact stands in marked contrast to the homogeneous system, which in the Tonks- 
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Girardeau limit had an infrared divergence p{k) = 0{k^^^'^). The seeming contradiction is owed 
to the fact that we deal with a bounded system, which cannot display true long-range order. 

Double well (h = 5) The momentum spectrum for a double well looks quite different from 
the start (g = 0.4): it exhibits two sidelobes. This can be explained by the symmetric orbital 
(f)oix) = c[(p{x — xo) + (p{x + xq)], which leads to a cosine-type modulation of p due to (f)o{k) = 
ccos{kxo)ip{k). These sidelobes get even more distinct as /i — oo (see Fig.lSt). With increasing 
repulsion (g = 4.7), there are two competing effects. On the one hand, the orbitals are flattened a 
little, which should result in a slightly sharper momentum distribution. It turns out, though, that 
the effect of fragmentation outperforms the former one even for tiny interactions: admixing an 
anti-symmetric orbital 0i adds a sin (/cxq) -type modulation, thus washing out the sidelobes as well 
as the central peak. Note that this effect is even more striking for h = 10, where it kicks in already 
for (jf = 0.2. In other words, the signature of the Gross-Pitaevskii regime in the harmonic trap — the 
initial sharpening of the k = peak — is lost in the case of a sufficiently pronounced double well. 

Along the lines of the remarks in the previous paragraph, we mention that the behavior for large 
correlations g is again universal as far as the k^^ tail for /c ^ oo is concerned. It also has a reduced 
peak for zero momentum, in accordance with the loss of long-range order found in Sec. lIV Al 

V. TWO-BODY CORRELATIONS AND DISCUSSION 

The focal point of our discussion so far has been the one-particle density matrix, and related 
quantities. However useful they are in studying fragmentation, they are by construction ignorant 
of an essential ingredient: the two-body correlations, which have of course been traced out in the 
definition of pi. Studying the (diagonal) two-particle density p2{xi,X2) may thus promise to yield 
an intriguing look behind the scenes of the one-particle picture. We will round off this section by 
commenting on the relation of our results to approximate methods such as multi-orbital mean-field 
theory and the two-mode model. 

A. Two-particle correlations 

For weak enough interactions, g ^ 0, the system is in an uncorrected state characterized 
by ^2(3^1,3:2) = p{xi)p{x2). The first effect of the two-body interaction is to distort the one- 
particle density governed by the Gross-Pitaevskii eq., while from some point on p2 will reflect 
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Figure 6: 2-particle density p2{xi,X2) for iV = 5 bosons in a double-well trap of barrier height h. Top row: 
h = 0, bottom row: h = 5; shown are the interaction strengths g = 0.4, 4.7, 10 from left to right. 

the correlations that are introduced by V{xi — Xj). As the pair-interaction energy is tr(Kp2) 
g J dxp2{x, x), keeping it low amounts to depleting the 'correlation diagonal' {xi = X2}. 

For g = 0.4 (Fig. a closer look reveals exactly that. In the case h = 0, a slight dip along 
the diagonal is formed as the density is being smeared out a little. In turn, for the double well 
(h = 5), the density is pumped into the off-diagonal peaks, which is indicative of a correlated 
state. Note that this 'dip' on the correlation diagonal is a feature of the two-body picture; in the 
one-particle density p = J dx2P2i-, X2) it is smoothed out and thus much less visible. As g is 
increased {g = 4.7), we are in the regime of fragmentation. Here the wave function develops 
marked minima at points of collision, Xj = Xj, as exemplified for two atoms in a harmonic trap 
IibI . This carries over to p2, where a characteristic correlation hole emerges, cutting the plot 
into halves. Despite that, the overall pattern still bares a resemblance to the non-interacting case. 
This changes when the system approaches fermionization, cf. g = 10. Here we encounter the 
simultaneous splitting into humps familiar from the one-particle density as presented in Ref. 1I20I1 . 
In this context, these wiggles signify that, if one boson resides at Xi, then any second one is likely 
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to be found at — 1 distinct spots X2, just not at X2 = Xi. The exact distribution depends on the 
trap, of course. For /i = 0, the checkerboard pattern is quite regular and has a larger amplitude 
about a; = 0, while for h = 5 the peaks are in a way packed into either well but suppressed in the 
center. 



B. Relation to approximation schemes 

Let us stop to wrap up what we have found and work out the key points by contrasting them with 
two well-known approximations. The transition from a macroscopic state 0f ^ to a fermionized 
state follows different pathways for different traps. For a harmonic trap, the lowest orbital retains 
its singular importance for any interaction strength. In a double well, as the barrier height tends to 
infinity, the route passes through a configuration fairly well approximated by two single-particle 
states (pa, i.e., a number state \N/2,N/2) (or generally a superposition of states \nQ,n[)). For 
higher g, in turn, the system will be fermionized, and the first single-particle orbitals will have 
dominant weights Ua. 

The stopover near the fragmented state \N/2, N/2) in the double well is the essence of the 



commonplace two-mode approximation (see, e.g., 112411 ). It can be recovered in MCTDH by re- 
stricting the number of single-particle functions to = 2 (see llllb . which sets the subspace to 
span{|no, n'l) | Xla'^a = Needless to say, it becomes only exact in the limit h ^ 00 and 
(7^0+. Never is it able to describe anything but that coarse-grained fragmentation into two sim- 
ple fragments, let alone the typical fermionization pattern in p(x), the correlation hole evidenced 
in p2{xi, ^2), or the short-range-driven tail evidenced in the momentum spectrum. 

The simple two-mode state above, |?2Q,n'^), is contained in a general multi-orbital mean-field 
theory, which yields the variationally best number state ll^ Issl . Among its most impressive 
successes has been a description of the evolution from the Gross-Pitaevskii regime to fermioniza- 
tion. The latter one was modeled by a state with = 1 (a < A^). We have pointed out [2(|] 
that this state emerges from our approach as a robust yet spurious solution if the Hilbert space 
is restricted to span{$j | < A^}. Also, we have argued that the number state gets all local 
quantities (the reduced densities as well as the energy) about right. By contrast, it is not designed 
to reproduce correlation- sensitive functions, such as the off-diagonal density matrix. Even though 
the use of several orbitals also serves to destroy off-diagonal long-range order (which hinges on 
the superposition of different orbitals, as delineated in lllBll . the incorrect populations {ria) make 
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Figure 7: The spurious momentum density p{k) as obtained within the restriction to n = N = 4 orbitals in 
the Tonks-Girardeau limit. Plotted are the results for the harmonic trap (h = 0) and the double well {h = 5) 
for large interactions g. The states correspond to number states |lo, . . . , lAf-i)- 

for afermionic rather than a bosonic momentum density, as displayed in Fig.|7] It resembles the 
spatial distribution and not the k = 0-peaked and long-range structure in Fig. |5l This reflects the 
fact that the mean-field approach cannot possibly recover the correct short-range behavior, which 
would require explicit correlations or, equivalently, the superposition of several single-particle 
configurations $j. Instead it only mimics the spatial separation. 

Lastly, one should be aware that for a state with Ua = I Va, the spectrum of pi is entirely 
degenerate; hence the eigenvectors (pa are only defined up to unitary transform (Uab)- In this light 
they might as well be thought of as spatially localized, as opposed to distorted oscillator functions 
(h = 0) or anti-/symmetric orbitals (h = 5). 



VI. CONCLUSIONS AND OUTLOOK 



The focus of this work have been correlation aspects of the numerically exact ground state 
of few atoms in different double- well traps. This way we have extended the more intuitive no- 
tion of the pathway from the weakly interacting regime via fragmentation to fermionization as 
described in our previous work lad. On the other hand, this article closes the gap between the 
non-interacting and the Tonks-Girardeau limit, for which aspects such as off-diagonal long-range 
order and the momentum spectrum have been understood, and shows how these very different 
cases connect. Our method is based on the Multi-Configuration Time-Dependent Hartree code, 
whose efficient variational approach allows us to compute the ground state to a high accuracy. 
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As one key result, we have highlighted the relation of the fragmentation process to the dimin- 
ishing of off-diagonal long-range order in the one-body density matrix. That mechanism has been 
explained in terms of the eigenvectors of the density matrix, the natural orbitals. These also allow 
us to relate the loss of long-range order to the momentum spectrum, whose markedly peaked struc- 
ture in the non-interacting case is stretched into a characteristic high-momentum tail for stronger 
interactions. Moreover, we show how the populations of all natural orbitals evolve, which not only 
illuminates how the fragmentation mechanism is altered as the trap is turned from a harmonic one 
to a double well, but also casts a light on the validity of the two-mode model. Finally, we have 
laid out in more detail the two-body nature of the correlations, which reflects in the formation of a 
'correlation hole' and culminates in the onset of the familiar checkerboard pattern. This goes well 
beyond the scope of approximation schemes inspired by a single-particle picture. 

With these investigations, the analysis of the ground state of few-boson systems in double-well 
traps may be considered complete. Future extensions of these studies appear obvious. For one 
thing, the addition of more wells to the trap is of interest. This touches on the question of the 
few-body analog of an optical lattice, and its related effects such as the superfluid/Mott-insulator 
transition. On the other hand, a physically thrilling situation would involve not only the ground 
state, but also excitations, and eventually as well looking into the dynamics of the system. Given 
the richness of these fields on the many-atom level, the detailed study of few atoms promises a 
wide range of applications. All these efforts may serve as a bridge toward a better control of 
ultracold few-body systems. 
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